Eco-evolutionary modelling of microbial syntrophy indicates the robustness of cross-feeding over cross-facilitation

Syntrophic cooperation among prokaryotes is ubiquitous and diverse. It relies on unilateral or mutual aid that may be both catalytic and metabolic in nature. Hypotheses of eukaryotic origins claim that mitochondrial endosymbiosis emerged from mutually beneficial syntrophy of archaeal and bacterial partners. However, there are no other examples of prokaryotic syntrophy leading to endosymbiosis. One potential reason is that when externalized products become public goods, they incite social conflict due to selfish mutants that may undermine any mutualistic interactions. To rigorously evaluate these arguments, here we construct a general mathematical framework of the ecology and evolution of different types of syntrophic partnerships. We do so both in a general microbial and in a eukaryogenetic context. Studying the case where partners cross-feed on each other’s self-inhibiting waste, we show that cooperative partnerships will eventually dominate over selfish mutants. By contrast, systems where producers actively secrete enzymes that cross-facilitate their partners’ resource consumption are not robust against cheaters over evolutionary time. We conclude that cross-facilitation is unlikely to provide an adequate syntrophic origin for endosymbiosis, but that cross-feeding mutualisms may indeed have played that role.


Examples of syntrophy
. Examples

Description of the interaction N/E Type Ref.
Cross-feeding between auxotrophic E. coli knock-out mutants lacking genes required for the biosynthesis of amino acids, nucleotides, and cofactors, as well as genes involved in glycolysis and respiration. In an experiment, 17% of random pairings out of 1035 combinations showed synergetic growth in co-culture of such mutants.
2 E CF flowthrough 1 Multi-member cross-feeding between engineered amino-acid auxotroph E. coli mutants. Biosynthetically costly amino acids tend to promote stronger cooperative interactions than cheaper ones. The evolved crossfeeding pair was also resistant against the invading autotroph wild-type strain. Besides pairwise, threemember communities (with double-auxotrophs) also showed synergistic growth when all three members were present. In another experiment, an initially 14-member systems undergo a drastic population shift toward a consortium dominated by four members, hence showing that microbes with multi-auxotrophic phenotypes can stably evolve, but only up to a limited size of the interaction network.
2 E CF flowthrough 5 Anaerobic methanotrophic archaea in a syntrophic relationship with sulfate-reducing bacteria. The partners rely on an efficient electron exchange via reducing equivalents or direct, tubular cell-to-cell connections. 2 N CF recycle 6,7 Successional colonization dynamics is often driven by cross-feeding where pioneer primary degraders enable late colonizers to feed on their by-products. For example, in particle-attached biofilms that degrade insoluble forms of organic matter concentrated on particles, community assembly is driven by a functional characteristic-based trophic structure: from narrow-range degraders to broad-range consumers, via resultant by-products of the degradation process. 2+ N CF flowthrough [8][9][10] The majority of bacterial species produce siderophores, a central component in the mechanism of iron scavenging. Siderophores are costly products secreted to make iron acquisition possible, yet the iron made available can be taken up by non-producers as well. After the iron is imported to the cell, the extracellular siderophore can be recycled (salvaged) or undergo hydrolysis. 11 Collective resistance can emerge in a community if an individual strain provides resistance for all members. It does so either by secreting extracellular enzymes to neutralize antibiotics (e.g. β-lactamase) or by neutralizing them intracellularly hence removing them from the environment. A particular form is the socalled cross-protection mutualism, in which sensitive and resistant strains in coculture can protect each other in a multidrug environment.

2+ N CX
2+ N CX [12][13][14][15] The yeast Saccharomyces cerevisiae secretes invertase enzyme to hydrolyze sucrose. Almost 99% of the end products, the monosaccharides, diffuse away benefiting other species nearby. 2+ N CX 16 The biofilm can provide protection against antibiotics originating from outside of the community, due to complex tolerance and resistance mechanisms, that may be more effective than resistance due to extracellular enzymes. The physical separation of the cells within the biofilm matrix from the external environment provides efficient protection against predation, invasion, bacteriophages, antibiotics, etc. [17][18][19] Surfactin bolsters surface spreading by reducing surface tension and promoting filamentous growth, benefiting all cells in the vicinity. It is costly to produce, hence acts like a common good, exploitable by cheaters.

2+ N CX
Pathogenicity in microbes is most often carried out by the secretion and acquisition of extracellular molecules that help to subvert the immune response of the host or facilitate the access to host nutrients, for example by triggering the outflow from the host. These products enable surrounding microbes to get access to energy rich nutrients of the host, acting as public goods that facilitate higher nutrient uptake for beneficiaries.
2+ N CX 21,22 Most bacterial species depend on cofactors acquired from the environment via cross-feeding, which they are incapable of producing due to their costly and complicated synthesis. It is predicted that the majority of bacteria harbor enzymes dependent on the cofactor corrinoid, which they have to acquire via cross-feeding, and which that enable pathways that significantly enhance their metabolic capacity.
2+ N CX 23 Rhodococcus species can degrade recalcitrant toxic pollutants in diverse environments, but their degrading efficiency is severely reduced by various environmental stress. Bacillus cereus exhibits strong stress resistance and can regulate pH. For example, Rhodococus ruber degrades a tetrahydrofuran providing acidic metabolites to B. cereus, which, in return, regulates pH and secretes nutrients essential for R. ruber.

N
hybrid, flowthrough 24 Up to five-membered cellulose-degrading community of various strains coexist stably under laboratory conditions exhibiting a complex, intertwined interaction network with mutualistic or commensalistic (syntrophy), inhibitory (as the accumulating intermediary product inhibits the catalytic activity of species), and competitive relationships. Both cross-feeding and cross-facilitative interactions are at play.

E
hybrid, flowthrough 25 Certain bacteria (e.g. Bacteroidales) live in the human intestine capable of breaking down polysaccharides extracellularly using glycoside hydrolases that in some cases are secreted outside of the cell. This breaking down is costly and results in diffusible public goods that benefit other species. For example, Bacteroides ovatus digests polysaccharides extracellularly, apparently benefiting only other species, and hence can be seen as an altruist act. The benefit of this act is enjoyed by nearby species, such as Bacteroides vulgatus, which reciprocates in various ways, including detoxification of inhibitory substances, or production of a growth promoting factor, initiating defense mechanisms and metabolism. This implies that this naturally evolved system involves both cross-feeding (by B. ovatus) and cross-facilitative interactions (by B. vulgatus). 2+ N hybrid, flowthrough [26][27][28][29][30][31][32][33][34][35] The bacteria Dehalococcoides mccartyi is the only known species that can completely degrade trichloroethene (TCE) from contaminated groundwater, but it depends on other species for a variety of exogenous compounds, such as hydrogen, acetate, corrinoids, biotin, and thiamine. Several other compounds, including hydrogen and CO, is hypothesized to take part in the syntrophic interactions. 2+ N hybrid, flowthrough 36 The consortium Chlorochromatium aggregatum consists of a motile central betaproteobacterium which transports its non-motile photosynthetic epibionts (Chlorobium) towards light while the epibionts provide energy for the host. Several other consortia are described with similar structures. These examples are the closest to syntrophic endosymbiosis and demonstrate the most sophisticated forms of pairwise prokaryotic ectosymbioses.

N
hybrid, flowthrough [37][38][39] In a metal working fluid-degrading community of low diversity (consisting of Agrobacterium tumefaciens, Comamonas testosteroni, Microbacterium saperdae, and Ochrobactrum anthropi), facilitative or neutral interactions dominate in a toxic environment. M. saperdae is highly dependent on by-products of other community members (cross-feeding), whereas cross-facilitative interactions arose by the other species removing toxic compounds initially present in the environment, enabling growth for the others. 4 N hybrid 40 Anammox communities mediate anaerobic ammonium oxidation (anammox) via complex interaction networks. Anammox granules consist of a mixture of cell aggregates and abiotic particles embedded within a matrix of organic extracellular polymeric substances (EPS). Within this medium, metabolite exchange takes place between heterotrophic and anammox bacteria. These include the products of anaerobic ammonium oxidation, nitrite loop as well as vitamins, and extracellularly degraded proteins. In particular, a Brocadia species demonstrating increased expression of genes involved in pathways for B1, B7 and B12 synthesis may indicate that this strain provides the above vitamins for the whole community.

Existence of a single, globally stable internal fixed point for cross-feeding without inhibition
We assume that the inhibitory effect of externalized waste product over the producer species (or anyone) can be neglected in the system described by (Eq. 4) in the main text. That is ℎ = 0 for all species (see Figure S1). The internal fixed points can be calculated analytically for two species (1 and 2), but their complexity prevents effective analysis of the system. Therefore, we rely on other means to ascertain that there is a single internal fixed point which always exists, and that it is stable. Figure S1. Model of cross-feeding without self-inhibition, with mutant species 3 . Arrows indicate consumption and production.
Squaring them indicates when one is guaranteed to be positive. In turn, the second ZNGI has the form 2 = − + √ 2 1 + 2 , which is qualitatively the same as the first one, just mirrored along the 1 = 2 identity line. That is, 1 = 0 and 2 = 0 defines convex and concave strictly monotonous curves in the 1 , 2 plane. Thus, they intersect each other at most in two points, of which at most one can be stable. The isoclines are depicted in Figure S2A. Figure  We now see that there is at most one internal fixed point. If this fixed point exists, then it is locally stable, due to the geometric arrangement of the ZNGIs ( 1 increases to the left of the blue curve and decreases to the right of it; conversely, 2 increases below the yellow curve and decreases above).

A B C
Assuming unrealistic parameter combinations, one can have two internal fixed points, of which only one would be stable (see Figure S2B). However, these cases imply that the per capita death rates are higher than the birth rates < (meaning that species are obligately dependent on each other), and/or an unrealistically high degradation rate of one of the byproducts.
Global stability can also be proved, using the Bendixson-Dulac criterion. This criterion states that, given some domain in the phase plane, the system 1 ′ = ( 1 , 2 ), 2 ′ = ( 1 , 2 ), and an arbitrary function ℎ( 1 , 2 ), if the divergence 1 (ℎ ) + 2 (ℎ ) is of the same sign everywhere in , then there can be no periodic orbits in that region. In this case, we can choose ℎ( 1 , 2 ) = 1/( 1 2 ) to immediately prove the impossibility of cycles for the whole positive quadrant: Since all parameters are positive, this expression is always negative in the whole positive quadrant. Thus, equating with this region, periodic solutions are impossible wherever 1 > 0 and 2 > 0.
Adding to this the fact that a two-dimensional differential equation system can never produce chaos, the only option left is that the locally stable fixed point must also be globally stable.
When does the fixed point exist? Nonexistence happens when the first isocline crosses the 2 = 0 line earlier than the second isocline (for an example, see Figure S2C). In the biological context, this means that species 2 has too great of a mortality to persist, and only species 1 prevails.

Numerical investigation of stable fixed points for cross-feeding with inhibition
We examine the fixed point distribution for a large number of random parameter combinations for the equation system (Eq. 4) in the main text. For each combination, we perturb the parameter values { 1 = 2 = 1, 1 = 2 = 0.1, 1 = 2 = 1, 1 = 2 = 1, 1 = 2 = 0.01, 1 = 2 = 0.1, 1 = 2 = 0.1, ℎ 1 = ℎ 2 = 0.5, 1 = 2 = 1} independently by adding a random value drawn from the normal distribution with mean 0 and standard deviation 0.5, ensuring that actual parameter values do not go below 10 -3 . We calculate the numerical Jacobian of the system with the given parameters, substituting in the numerically calculated fixed points as the densities. For each fixed point, we check whether the eigen values of the Jacobian are all negatives. After extensive numerical testing (more than 100 000 independent random parameter combinations), we have not found a parameter combination which has more than 1 internal stable fixed point.

Existence of a single, globally stable internal fixed point for cross-feeding with inhibition
We now include self-inhibition of waste materials on species into the dynamics (see Figure S3). We follow the same analysis as before, in SI 3. The per-capita growth rates of the subsystem of species 1 and 2 are: At the internal fixed point, where 1 , 2 > 0: 1 = 2 = 0. Due to the complexities of the dynamics with self-inhibition, we cannot analytically or geometrically prove the existence of a single globally stable internal fixed point. However, based on the extensive numerical investigations we have performed (see above), it is strongly believed that there is at most one internal stable fixed point of the 2-species system.
The same divergence formula can be used as before, to prove the impossibility of cycles for the whole positive quadrant ( 1 , 2 > 0), according to the Bendixson-Dulac criterion, assuming 1 ′ = ( 1 , 2 ), 2 ′ = ( 1 , 2 ): Since all parameters are positive, this expression is always negative in the whole positive quadrant. Thus, periodic solutions are impossible wherever 1 , 2 > 0. Therefore, the locally stable fixed point must also be globally stable.
The zero net growth isoclines behave similarly to the system without self-inhibition, depicted in Figure S4. Again, biologically unrealistic parameters may result in the disappearance of the internal fixed point or the appearance of a second, instable internal fixed point, like without inhibition.

Condition of invasion for a third species for cross-feeding
We would like to know whether a mutant species can invade a presumably ecologically stable pair of cross-feeding species in case there is no inhibition by the products . Therefore, we analyse the invasion dynamics of the model. We assume that species 1 and species 2 are present initially (Eq. S1), and we ignore self-inhibition of metabolites (ℎ = 0, ∀ ).
Above we have proven that there is only one stable positive fixed point of this system. Let us denote the fixed points of species 1 and 2 as ̂1 and ̂2, respectively. Because of the definition of fixed points, it is true that: Let us use the following notations: 3 = 2 + ∆ , 3 = 2 + ∆ , 3 = 2 + ∆ , where the differences between the parameters , and (∆ , ∆ and ∆ respectively), can be negative positive or zero. Substituting (Eq. S3) into the dynamical equation of 3 and taking the linear approximation of it by assuming that 3 ≪ 1: In case of ∆ = ∆ = 0, and ∆ > 0 that is when withholding 2 increases death rate of species 3 comparing to species 2, then species 2 can invade in the species 1, 3 subsystem. That is, species 2 and 3 cannot coexist.
Another interesting case is when ∆ , ∆ > 0 at the same time. In this case, not secreting the byproduct 2 means extra cost for the mutant species 3. However, it utilizes the byproduct 1 of species 1 more effectively than species 2. According to (Eq. S4), species 3 can spread if: and similarly, species 2 spreads in the species 1, 3 subsystem if: It is easy to show that ̂1 * < ̂1and assume that ̂3 ≥̂2 (which, depending on the parameters can be valid or not), then the relations (Eq. S6) and (Eq. S7) can be satisfied simultaneously. In this case, species 2 and 3 mutually invade each other, thus all three species will be in coexistence. If only (Eq. S6) is true and (Eq. S7) is not, then the cheater mutant excludes the cooperator one. However, the evolutionary history depends on the trade-off between and . It is natural to assume that an increased efficiency in consumption entails an increase in death rate, so these two variables are in positive trade-off. This is a reasonable assumption, as the increased performance of any metabolic machinery of incurs an energetic cost that either decreases reproduction or growth rate or increases the death rate. The evolution of these correlated traits via adaptive dynamics is discussed in SI 6 and SI 8.
In turn, we examine the invasion dynamics of a third species when there is self-inhibition of products (ℎ > 0, ∀ ). Since there is only one positive stable fixed point in the model of cross-feeding with selfinhibition (SI 5), we can use the same analysis as we did for the case without inhibition. Again, we consider the species 1, 2 subsystem. From the definition of fixed point, it follows that where ̂1 and ̂2 are the positive fixed point of the dynamics. Using the similar notations as above and denote ℎ 3 = ℎ 2 + ∆ℎ, we receive that the in the linear limit the mutant species 3 dynamics is: 3 .

Adaptive dynamics for cross-feeding without inhibition
We assume the following trade-off between the efficiency of consumption of the byproduct ( ) and mortality ( ): where is the steepness of the sigmoid transition, ̅ is the expected mean of the trait value , and is a scaling factor for positive values of .
We assume that the underlying trait, , can have any value, even negative, but as increases, both dependent traits and increase. For a visual explanation of the trade-off, see Figure S5. The modified rate equations are: where the growth of species 1 (the unchanging resident) is governed by the first equation, while the second equation describes the growth of any mutant of species 2 (itself included), and is the total number of species in the system.
First, we examine the system when species 1 is not allowed to mutate, only species 2 (see Figure S6). Figure  Starting each simulation from a different initial trait value for species 2 shows the different evolutionary trajectories that lead to one of two equilibrium trait values (see Figure S7). We use the per capita growth rates of species 1 and 2 (Eq. S2) to check if at the equilibrium, no mutant of species 2 can invade the system. The invasion growth rate 2 * of a mutant of species 2 with trait 2 * is:

A B
At equilibrium, d d = 0 and d d = 0 for both species, with equilibrium value =̂. Analysis reveals that the growth rate is always negative for a mutant when species 1 and 2 are in equilibrium (see Figure S8). Note that species 1 is not stable evolutionary in this case but is kept fixed. In case both species can evolve, we define two general growth equations for the two mutant classes with respective abundances 1 and 2 : where is the total number of species belonging to the th mutant class. The two class of species, when their initial trait values allow, evolve towards better cross-feeders (see Figure S9). Figure S9. = 0.01, = 10 −2 , = 10 −3 , = 1000, 1 = 0.2, 2 = 0.3, 1 = 1, 2 = 0.1, 1 = 1, 2 = 0.1, 1 = 1, 2 = 0.1, 1 = 1, 2 = 1, 1 = 2 = 0.1, 1 = 1, 2 = 1, ̅ = 0.5, = 0.2, = 0.1}. While species of have higher birth rates than species of , they are present in less amounts at equilibrium (opacity of points correspond to relative equilibrium density). Nevertheless, they can evolve to better trait values than species of .

Adaptive dynamics for cross-feeding with inhibition
We take species 1 to be a static resident and only allow species 2 to mutate: where is the total number of species. In case both species can evolve, we define two general growth equations for the two mutant classes with respective abundances 1 and 2 :  where is the total number of species belonging to the th mutant class.
Simulations where only one or both species can evolve, starting from different initial mutant trait values are shown in Figure 4 in the main text.

Coexistence of species 1 and 2 for cross-facilitation
We study the case when only species 1 and species 2 are present (see Figure S10). The dynamical equations of the concentrations can be estimated below and above by: ) for all positive initial densities. This means that species 1 and species 2 coexist. Naturally, we have no information about the number of fixed points and their characteristics by this analysis.

Numerical investigation of stable fixed points for cross-facilitation
We examine the fixed point distribution for a large number of random parameter combinations for the equation system (Eq. 6) in the main text. For each combination, we perturb the parameter values { 1 = 2 = 1, 1 = 2 = 0.1, 1 = 2 = 1, 1 = 2 = 1, 1 = 2 = 0.01, 1 = 2 = 0.1, 1 = 2 = 1, 11 = 12 = 21 = 22 = 0.1} independently by adding to each parameter a random value drawn from the normal distribution with mean 0 and standard deviation 0.5, ensuring that actual parameter values do not go below 10 -3 . We calculate the numerical Jacobian of the system with the given parameters, substituting in the numerically calculated fixed points as the densities. For each fixed point, we check whether the eigen values of the Jacobian are all negatives. After extensive numerical testing (more than 100 000 independent random parameter combinations), we have not found any parameter combination which has more than 1 internal stable fixed point.

Existence of a single, globally stable internal fixed point for cross-facilitation
The per-capita growth rates of the subsystem of species 1 and 2 are: ( 1 1 + 1 1 ) 2 2 − 2 1 1 1 11 1 1 ( 1 1 + 1 1 ) 2 2 − − 1 1 2 11 1 2 ( 1 1 + 1 1 ) 2 2 − which may be negative, depending on a complex combination that cannot be simplified (assuming 1 , 2 > 0 and all parameters are positive). While we did not manage to prove that only a single internal stable fixed point exists in this system, we have never found any other dynamical outcomes despite extensive numerical investigations (SI 10). We therefore conclude that even if multiple stable states are a possibility, they are highly unlikely to be realized and one can proceed as if the internal equilibrium point was unique. A characteristic phase plot is shown in Figure S11.

Adaptive dynamics for cross-facilitation
We use the same trade-off functions as before (Eq. S8), now between the traits mortality ( ) and production rate ( ), as we assume that production is costly, and the higher the production rate, the higher the death rate is (see Figure S12 Assuming that both species can evolve, we define two general growth equations for the two mutant classes with respective abundances 1 and 2 , modifying the dynamical equations of (Eq. 6) in the main text: where is the total number of species belonging to the th mutant class.
Results for adaptive dynamics when both species can mutate are shown in Figure 5 in the main text.
If is close to zero (but not identical), several species can live together without exclusion (see Figure S13). This is because species with ≤ 0, their mortality is effectively zero ( ≈ 0), and marginally small (having no effect), which leads to a degenerate case where species can coexist in a neutral equilibrium. Figure

Adaptive dynamics for another model of cross-facilitation
In case of cross-facilitation of Figure 1B (and (Eq. 5) in the main text), products 1 and 2 have identical effects: they both improve the consumption of both resources 1 and 2 . Accordingly, 1 and 2 can actually represent the same molecule, instead of two different products. Here, we modify the cross-facilitation case of (Eq. 5) to model the situation where product improves the consumption of only, but makes the resource available to any species. In this setup, we effectively differentiate between 1 and 2 , as they cannot be the same molecule (see Figure S14). Ignoring species 3 for the cross-facilitation case, the dynamical equations are: Parameters are different for the two mutant classes { = 0.01, = 5 • 10 5 , = 10 −2 , = 10 −3 , = 450, 1 = 2 = 1, 1 = 2 = 0.1, 1 = 1, 2 = 0.1, 1 = 1, 2 = 0.1, 1 = 2 = 1, 11  14. Adaptive dynamics for a hybrid cross-feeding, cross-facilitation system Next, we investigate the case where one of the species is a cross-feeder and the other is a crossfacilitator, as depicted in Figure 7.4 in the main text and detailed in Figure S16. Figure S16. A possible hybrid case of syntrophy in which one species is a cross-feeder, the other is a crossfacilitator. Species 3 is a mutant of the cross-facilitating species 2 that does not produce anything but benefits from the catalytic effect of 2 (dashed arrows).
Accordingly, species 1 enjoys a benefit of the enzyme of species 2, but plus the benefit of species 2 consuming its waste. On the other hand, species 2 enjoys the benefit of eating the waste of species 1 plus the benefit of its own external enzyme. The dynamical equations are as follows: In turn, the metabolite dynamics read: We assume fast dynamics for the by setting Substituting fast resource dynamics (of (Eq. 1) in the main text) and fast metabolite dynamics back into the dynamical equations yields: d 1 d = 1 ( 1 ( 1 − 1 1 ) − 1 − ℎ 1 1 1 1 + 2 2 + ( 1 − 1 1 ) 12  ).
Evolution mutates a unidimensional trait in case of cross-facilitators. Their mortality and enzyme production rate depend directly on through the same trade-off as before, as defined in (Eq. S9). However, in this setup, waste-consumption efficiency 2 is also a property of cross-facilitators, and not of cross-feeders. We keep it as a static property. As a minor consequence, cross-feeders cannot consume the waste of other species hence they only have a single property that could depend on a possible evolutionary trait for the cross-feeder class. For sake of simplicity, the cross-feeding species ).